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Abstract 

Many interfacial phenomena in physical and biological systems are dominated by high order geometric 
quantities such as curvature. Here a semi-implicit method is combined with a level set jet scheme to handle 
stiff nonlinear advection problems. The new method offers an improvement over the semi-implicit gradient 
augmented level set method previously introduced by requiring only one smoothing step when updating 
the level set jet function while still preserving the underlying methods higher accuracy. Sample results 
demonstrate that accuracy is not sacrihced while strict time step restrictions can be avoided. 
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1. Introduction and Overview 


Level set methods are one of the most popular interface capturing approaches with a wide variety of 
applications ranging from fluid flow [l| , image segmentation Q , materials science and structural topology 
optimization P). In the level set approach an interface r(t) is represented as a zero level set of a smooth 
function t): 

T{t) = {x:(j){x,t) = 0}. (I) 


The evolution of the interface is implicitly tracked by evolving 4’{x, t) due to an underlying velocity field u, 
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Level set methods are robust in handling topological changes and also enable easy calculation of geometric 
quantities of the interface. Eor instance, the unit normal pointing outward from the interface and the 
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curvature of the interface can be expressed in terms of 4>{x,t): 
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The accuracy of the level set method can be improved by including the gradients of the level set and 
advecting them in a fully coupled manner. This method is known as the gradient augmented level set and 
was first introduced by Nave et al. as a globally third order accurate method that uses p-cubic Hermite 
interpolation with an optimally local stencil Q . Benefits of incorporating gradient information also include: 
1) representation of subgrid structures, 2) simple and accurate approximation of the curvature and 3) simple 
implementation of adaptive mesh refinement with optimal locality Q] . When us^ for multiphase flow these 
benefits result in increased volume conservation over the base level set method 


0 . 


Jet schemes, developed by Seibold et ah, are a natural extension of the original gradient augmented 
level set method and allow for the creation of schemes with an arbitrary high order of accuracy [8|. The 
higher order of accuracy is accomplished by tracking a suitable jet of the solution, i.e. appropriate derivative 
information in addition to function values. The higher order spatial derivatives are utilized to construct 
high order Hermite interpolations. Jet schemes of up to fifth order accuracy were developed and shown to 


obtain comparable results to that of WEN05 [sj. Chidyagwai et al. demonstrated that jet schemes are 
computationally more efficient than comparable WENO schemes while possessing the optimal locality of 


discontinuous Galerkin methods 


0 . 


Many problems concerning moving interfaces require computation of higher order spatial derivatives 
of the level set. For example, the surface tension driven flows of immiscible fluids contain nonlinear and 
nonlocal higher order terms due to the Laplace-Young condition at the interface 0 1^, 111 12, 13|. As another 
example, in modeling the hydrodynamics of lipid vesicles the bending resistance exerted on the surrounding 
fluid is proportional to the second derivative of the curvature [l^. This involves the computation of up 
to the fourth order derivative of the level set function which makes the evolution of the interface front 
extremely stiff. The gradient augmented level set method and jet schemes were originally developed for 
linear advection equations and are not able to maintain a smooth level set field without strict time step 


restrictions 


Q- 


The semi-implicit formulation of an Eulerian level set method was first introduced by Smereka and allows 


for time steps on the order of the grid spacing [l^. This method is based on extracting the linear portion 
of the advection equation and treating that implicitly while using an explicit form of the nonlinear part. 
Semi-implicit methods are a powerful alternative to explicit integration methods which suffer from strict 
time step conditions and fully implicit methods which require solving a non-linear system of equations for 

every time step HQ. 
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A combination of the standard semi-implicit method and the original gradient augmented level set method 


Q- 


was recently reported under the name of the Semi-implicit Gradient Augmented Level Set (SIGALS) 

The method provides more accuracy through tracking the level set gradients and alleviates some of the 
nonlocal behavior of the original semi-implicit level set method. Another closely related formulation of 
this method has been successfully implemented in the simulation of vesicle electrohydrodynamics and the 
numerical experiments show the robustness of the semi-implicit formulation in capturing the fast and long- 
range deformation of the vesicle in strong DG electric fields [^, 1^. In another recent work a semi-implicit 


method was introduced to treat the numerical instabilities of partial differential equations in general 


li- 


In this work the level set jet scheme is extended to model interface evolution due to numerically stiff 
velocity fields. Unlike the SIGALS method, which requires four smoothing operations per time step for a 
three-dimensional problem and updates the gradients through an analytical differentiation method, here a 
different approach is introduced to keep the level set and gradients smooth with less computational cost. 
Additionally, the derivatives are updated through a simpler and more effective e-finite difference method. 
Gompared to the original SIGALS method, which uses an analytic updating method for the derivative 
fields, this e-finite difference method extends more easily to higher level-set derivative fields. It has also 
been previously demonstrated that the accuracy of the e-finite difference method and analytic method is 
identical, while the e-finite difference method has a lower computational cost Q. 

The remainder of this paper is organized as follows. In Section [5] the numerical background information 
is presented. This includes the base level set jet scheme and the original level set semi-implicit formulation. 
In Section |3] the new semi-implicit jet scheme (SemiJet) is described. Gonvergence and stability analysis 
using sample two and three dimensional results are provided in Section |4l A brief conclusion and future 
work is given in Section [SJ 


2. Numerical Background 

The method presented in this paper is built upon the original level set jet scheme and a semi-implicit 
formulation for Hamilton-Jacobi equations. In this section some background information on the original 
level set jet scheme and the semi-implicit level set method are given. A brief mention of velocity extension 
and reinitialization is also presented. 

2.1. Level Set Jet Schemes 

Gonsider the advection of an interface using the level set method on a fixed grid. The level set function is 
only known at grid points. To obtain interface information away from grid points some type of interpolation 
must be used. The idea behind level set jet schemes is to advect some or all of the information needed to 
obtain high-order interpolation functions in a local manner without the need to use wide stencils Q. A 
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grouping of this information, called a jet, consists of the level set function and certain derivatives of the level 
set. For example, using a jet which consists of the level set function, (/), in addition to gradient vector field, 
(j)x and <j)y, and the first cross-derivative, (j^xyi it would be possible to define a cubic Hermite interpolant on 
a two-dimensional Cartesian grid without the need for derivative approximations. 

A critical component of any level set jet scheme is the coupled advection of the jet information. The origi¬ 
nal gradient-augmented level set method and the more recent level set jet schemes both use a semi-Lagrangian 
Mproach: they evolve values on a fixed Eulerian grid by tracing characteristics using a Lagrangian method 
a y. Consider first the advection of the level set function alone. In this case the time-evolution is given by 
Eq. (12). Using the material derivative form of the evolution equation it can be written to first order that 

D4> (j)n+l - (j)'^ 


= 0 , 


(5) 


Dt At 

where (pn+i is the level set value at a grid point x at time tn+i, (t>n is the level set value at the departure 
(or foot) location, and At is the time step. The departure location can be obtained from 

Dx{t) 


Dt 


= u{x[t),t). 


To first order this can be written as 


x‘^ = X — At u{x, tn). 


The updated level set value can then be obtained by evaluating an interpolation function P^, 

0n+l — — P(f){x ,tn)- 


( 6 ) 


(7) 


( 8 ) 


Higher order method can be developed and will be discussed later. 

A critical component of Jet schemes is the advection of the additional information. This must be done 
such that derivatives of the level set do not de-couple from the base level set function. There are two 
basic methods to accomplish this. The first is called analytic differentiation and was used in the gradient- 
augmented level set method Q. In the analytic differentiation method spatial derivatives of the level set 
evolution equation are used to evolve of the derivative fields. For example, to advect the gradient of the 
level set, il’ = an advection equation can be obtained by taking the gradient of Eq. (ED, 


Dcj) 

15t 


— Vu ■ ■0. 


(9) 


The analytic differentiation method does not work well with higher-order derivatives as the resulting 
evolution equations can become quite cumbersome. An alternative is the e-finite difference method. At 
every grid point a local Cartesian sub-grid is constructed with a spacing of e <?; h, where h is the grid 
spacing. This defines a grid given by x^ = x + qe for q G { —1,1}^ where p is the dimension of the 
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Figure 1: The subgrid used in the e-finite difference update method for a Jet. 


simulation, either p = 2 or p = 3, Fig. [TJ At each subgrid point an updated level set value, (jfl^ is obtained 
using Eqs. © and ([H|). 

The updated level set values on the subgrid can be used to update the derivatives at grid points by 
using finite difference approximations. In two dimensions values up to (pxy can be obtained via the following 


finite-difference approximations: 

4> -I- -I- , (10) 

(jjx , ( 11 ) 

(t^y , ( 12 ) 

=4^ ~ . (13) 


The overall accuracy of the scheme is O(e^) -I- 0{d) for Eq. (ITO . O(e^) -1- 0{5/e) for Eqs. (ITTl) and (IT^ . 
and 0{e^) -I- 0{S/e^) for Eq. (HTTI) . where (5 is machine precision [8|. An optimal choice is thus e = 
which yields an overall accuracy of 0(5^/^) in all cases. Note that this particular choice for e ensures that 
the overall errors are dominated by other terms, such as the time discretization, and not the calculation 
of derivatives fields during updating. While, in general, the derivatives in a level-set description of an 
interface are discontinuous in a finite number of points for a discretized domain, these discontinuities are 
typically away from the interface and situations where they occur near the interface are temporary for 
moving interfaces. These discontinuities have not been observed to be detrimental to the overall accuracy 
of Jet schemes in the past @,01. 

The major advantage of the e-finite difference method over the analytic differentiation method is the 
applicability to higher-order derivatives. Higher-order derivative values can be calculated by expanding the 
subgrid stencil and using appropriate finite difference approximations. The small subgrid size ensures that 
results are accurate. Additionally, Chidyagwai et al. showed that e-finite difference method are faster than 
analytical differentiation and that the method is more efficient than WENO while possessing the accuracy 
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of discontinuous Galerkin schemes 


2.2. The Semi-Implicit Level Set Method 

One issue with all interface tracking techniques is their instability if the underlying velocity field is 
numerically stiff. A classic example is mean curvature flow where the velocity of the interface is given 
by M = —KTi. From Eq. 0 it is clear that the curvature is a second-order derivative of the level set 
function, which would require a time step restriction of At = 0{h?) for numerical stability if an explicit 
time discretization is used. 

The semi-implicit level set method was developed to circumvent this time step restriction [l^ . If the 
evolution is relatively simple, such as mean-curvature flow, then it is possible to extract the linear portion of 
the level set evolution equation. This linear portion would be done implicitly while the remaining non-linear 
portion is done explicitly: 

rh 1 -1 — rh _ \/rh 

V(||V</.„||). (14) 


4’n-\-l <Pn r-r2 1 ^ 4^n 

= V (pn+1 - 


At ' 

In cases where it is not feasible to explicitly extract the linear portion a smoothing operator can be added 
to the evolution equation Ql Q,[3, 

4^n+l 4^‘ 


At 


T '^n ' ^4^n -t- Z/((()^) — 0, 


(15) 


where the L((j)) is a linear damping operator which can be chosen based on the stiffness of the velocity held 
u. If the order of the velocity held is even, a common choice for the damping operator is 
where m matches the order of the velocity held and /3 is a constant. Note that this stabilization is added 
after normalization and therefore P is dimensionless. If the system is dimensional, then this /3 is akin to a 
diffusivity parameter. If the velocity order is odd, then a more appropriate choice is the Hilbert transform 
of the m*^-order spatial derivative, T(d)) = /3TL{V^(f)). This transform ensures that the right scaling is 
observed in the frequency domain [19|. 

An explanation of why Eq. (1151) aids in stability has been provided in Ref. [l^ and is briehy explained 
here. Consider the simpler problem of Ut = au^x on a one-dimensional domain with periodic boundary 
conditions. As the system is of second-order, an appropriate linear damping operator has the form of 
= Puxx, so the modihed evolution equation is now Ut — jSuxx + Puxx = ctUxx- Consider the perturbation 
of a single Eourier mode, u{xj,tn) = where ^ is the amplification factor, h is the grid spacing, At 

is the time step, and k is the Fourier mode. Using the perturbation in a first-order in time and second-order 
in space discretization and simplifying results in the following expression for the amplification factor. 


^^‘ = 1 


aAt ( ikh 

~w 


^ — ikh 


-2) 


1 _ ^ (^e^kh 


:}—ikh 


- 2 ) 


Using the relation -|- e ^kh _ 2 — —4sin^(/ifc/2), the amplification factor is now 


4aAt 


= 1 - 


^sm^{hk/2) 
l + ^sin2(hfc7^' 


(16) 


(17) 
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A stable scheme is one where < 1, which corresponds to the relation 

,^ At , 2 f kh\ 

(4a - 8/3) p- sin j < 2. 

As the sin^ term has a maximum value of 1, stability is ensured when 

(4a - 8/3) ^ < 2. 


( 18 ) 


(19) 


When /3 = 0, which corresponds to no stabilization, the standard time-step restriction for second-order 
differential equations is obtained, aAt/h? < 1/2. Any value of /3 > 0 stabilizes the scheme as the condition 
now becomes aAt/h? < 1/2 -|- 2/3At/h^. In this specific case, when /3 = a/2, the system becomes is 
equivalent to a Crank-Nicolson discretization and is unconditionally stable. For general advection, the value 
of /3 = 1/2 has been shown to work well in most cases [l^. Additionally, it has been shown that any scheme 
can be stabilized, even if the linear damping operator is of lower order than the velocity field [l^ . 

Finally, a word should be said about the boundary conditions for the semi-implicit update. If the domain 
is periodic then no additional boundary conditions are required. For non-periodic boundaries a logical choice 


is to use homogeneous Neumann boundary conditions 


2C 


it would be necessary to use an interface capturing technique 


211. In the case of complex or curved boundaries 


2al23| . Note that if level set reinitialization 


is used, as described below, the particular choice of boundary conditions can play little role in the evolution, 
as the level set function far away from the interface is replaced with a different, but equivalent one. 


2.3. Level Set Reinitialization and Extension 

While any level set function where a contour describes the interface can be chosen, a typical choice is 
to use a signed-distance function where the interior is given by (/ < 0 while the exterior has (/ > 0 and 
|(/)(a:)| indicates the distance from x to the interface. Even if a signed-distance function is chosen as the 
initial level set it can not be guaranteed that it will remain so over time. Typically level set reinitialization 
is performed periodically to return the level set back to a signed distance function. One possible method is 
to use a PDE-based reinitialization technique 


d(j) 


-sgn((/o) (||V(/)|| - 1) = 0, 


( 20 ) 


where r is a pseudo-time and sgn((/o) is the sign function of the original level set. Another possible method 
is to reinitialize by explicitly calculating the closest point from a point in the vicinity of the interface to the 
interface . While slower computationally, this explicit calculation of the closest point has the advantage 
of higher accuracy than the PDE-based method and will have additional benefits for Jet schemes, as will be 
explained later. 


In many situations the velocity is only defined at the interface. To advance a level set function the velocity 
needs to be determined, at a minimum, in a region surrounding the interface. An extension algorithm can 
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accomplish this by extending a function in the direction normal to the interface. One method to do this is 
to solve a hyperbolic PDE in the form of 


—+sgn(^)n ■ Vg = 0, (21) 

where q is the quantity to be extended from the interface. This PDE can be solved for a few iterations 
using a pseudo-time step of Ar « h/2 to extend the quantity g in a region near the interface. An alternative 
to the PDE-based extension algorithm is to explicitly extend a quantity by interpolating the quantity q at 
the closest point which could have been determined during a reinitialization procedure. In practice both 
methods work well. 


3. The Semi-Implicit Level Set Jet Scheme (SemiJet) 

The idea behind the SemiJet method is to combine the increased accuracy of the standard level set jet 
method with the stability of the semi-implicit level set method. The general outline has three steps: 1) a 
standard semi-Lagrangian, semi-implicit update, 2) determining the effect of smoothing, and 3) propagation 
of the smoothing effect to derivative updates, see Fig. [2] for a graphical representation. Without loss of 
generality, only even-ordered velocity field will be considered. Therefore, the linear damping operator will 
have the form of If another linear damping operator is needed, such as when the velocity 

field is dominated by odd-order derivatives, then is simply replaced by an appropriate operator. 



(a) (b) (c) 


Figure 2: The three basic steps of the SemiJet method, (a) A standard semi-Lagrangian semi-implicit step which updates grid 
points, (b) The effect of smoothing is applied to the sub-grid locations by way of a smoothing source term, (c) The subgrid 
points are updated using a semi-Lagrangian advection with the smoothing source term, allowing for the remaining Jet fields 
to be updated. 


The first step is a standard semi-Lagrangian semi-implicit step on the base level set function cj). Begin 
by writing the semi-implicit level set equation in semi-Lagrangian form, 

-/3V™<^„+i+/3V™0 = O, 


D(t) 

iDt 
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where (j) is an approximation to (pn+i- A first-order discretization in time results in 


(t>n+l - (t>i 

At 


/3V™(/.„+i+/3V™0„ = O, 


with (f> = (pn- This can be solved through a multi-step process given by 


( 23 ) 


x‘^ = X — At u{x,tn), 


(24) 



At 


{x‘^,tn) , 


(25) 

(26) 


where (/)(( is the tentative (departure) level set value and {x'^,tn) is the interpolant of the level set value 
at time evaluated at the departure location x'^. 

A second-order method can be obtained by writing 

~2At ^ = 0’ (27) 

where p = 2pn — Pn-i- The departure value of the level set at time is calculated by 


'^n+l — ‘2Un (^S) 

xi = X — At iin+i, (29) 

X^ = X - ^At{Un+l + Pu(xi,tn)) , (30) 

4‘n ~ Pip I (31) 


where Pu is the interpolation function of the velocity field. The departure value at time tn-i can then be 
obtained with 


= x- 2AtP.a , 

(32) 

= Pp ■ 

(33) 


A smooth level set held can now be obtained by solving Eq. (E3. Note that, while it is possible to use 
higher order time discretizations, the results presented below use second-order accurate approximations to 
the smoothing held, V™. In many cases the velocity held is normalized to be 0(1), and it is typically 
not advisable to have the interface move more than one grid spacing per time step. Therefore, a situation 
where At Ri h, where h is the grid spacing is a standard situation and thus it is not advantageous to use a 
higher-order time discretization as the method would be limited by the spatial discretization. 

At this point a smooth and updated level set held has been determined but the remaining derivatives 
helds in the Jet have not yet been updated. The application of the smoothing step with /3 > 0 has the effect 
of modifying the evolution of p, with the magnitude of this deviation from the /3 = 0 (e.g. no smoothing) 
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case given by /3V™^„+i — j3V^(j). Conceptually, if this deviation was known at the beginning of a time step, 
there would be no need to perform a semi-implicit update, as the corrections needed to maintain a smooth 
interface would already be known. 

This concept is used to update the remaining fields in the Jet. Capture the effect of smoothing by 
defining a smoothing source term, 

(34) 

where both (pn+i and 4> have already been determined. The method now proceeds in a similar manner to 
the standard level set Jet scheme, with the only modification that the advection now has a source term. At 
each grid point construct a sub-grid of spacing e <S h: x'^ = x + qe ior q € { — 1,1}^, see Fig. |2(b)[ At each 
sub-grid location calculate an updated level set value, For example, a first-order scheme would use the 
following steps: 


x‘^ = x‘^ - At Pu{x‘^,tn), 


(35) 



At 


{x‘^,t„) , 

Ps ix'i ), 


(36) 

(37) 


where Ps is the interpolant of the smoothing source term given in Eq. (I34p . The second-order scheme 
proceeds similarly, 


u (®^; ^n) Pu (^*^; l) 5 

Xi=x'^ - At iln+l, 

x^ = x'‘ - ^At {Un+1 + Pu {Xl,tn)) , 

a;^_i = x'^ - 2AtPu {xn,tn ), 

4’n—l ~ P<t> {x^_i,tn — l) , 

3^q - 4</.(( + 


2At 


= Ps {x"^) ■ 


(38) 

(39) 

(40) 

(41) 

(42) 

(43) 


In both the first- and second-order cases the smooth source term is evaluated at the sub-grid point x'^. 
Once the sub-grid level set values have been determined the updated derivative fields at grid points can be 
calculated using finite difference approximations such as Eqs. (HH)-®, or similar. Note that in the SemiJet 
method the level set function is not given by Eq. (nni, but is instead the result of the smoothing step, either 
Eq. (HH) or (l?7l) . 

A summary of the SemiJet method can be written as: 


1. Determine a tentative updated level set field using a standard semi-Lagrangian level set method. 

2. Solve the linear system arising from a seini-implicit update. This produces the updated level set field. 
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3. Calculate the effect of smoothing by way of a smoothing source term. 

4. Use the smoothing source term to determine updated level set values on the sub-grid. 

5. Update the remaining fields in the Jet using finite difference approximations on the sub-grid. 

From this summary is it clear that only one linear system needs to be solved per time step. This is an 
advantage over the original SIGALS method which required a linear system for each component of the Jet. 


4. Numerical Results 


In this section sample two- and three-dimensional results using the method described in Sec. [3] are 
presented. Unless otherwise specified a constant of /? = 0.5 is chosen and second-order accurate time 
integration is used. In this work only curvature flows are considered and thus the smoothing operator is 
set to the Laplacian, m = 2. The sub-grid spacing was set to e = 10“^ while periodic boundary conditions 
were used in all cases. The smoothing operator is discretized using second-order accurate isotropic finite 
differences 


26 


27| . The linear system of equations resulting during updating the level set function are 
solved using Flexible Generalized Minimal Residual (FGMRES) method with Geometric Algebraic Multi- 
Grid (GAMG) as the preconditioner, both provided by the PETSc library 


28 


29|, with a relative 


tolerance of the solver set to 10“®. 

All required interpolation functions are cubic interpolants. The interpolation coefficients for the velocity 
field, Pu, and smoothing source term, Pg, are determined using a standard 16-node stencil in two-dimensions 
and 64-node stencil in three-dimensions. The stencil for the interpolation coefficients for the Jet, P^, are 
determined based on the order of the Jet. The Jet structures here will be referred to as either a 0-Jet 
or Partial 1-Jet (Pl-Jet). These names refer to the number of derivative fields which are available for 
interpolation functions. A 0-Jet consists of just the level set field, (j). Thus a 0-Jet will require either 16- 
or 64-nodes, depending on the dimension of the simulation, to determine the cubic interpolant. A Pl-Jet 
includes the gradient field, ■0 = {(pxT 'Py^ 'Pz)- To compute the Jet interpolation function the first cross¬ 
derivative are also required. For a Pl-Jet these are computed using the cell-based approximation Q. In 
two-dimensions the values of pxy are first calculated at edge-centers. These are then interpolated to the 
cell corners (the grid points) to obtain the necessary value for the interpolant. In this way only information 
from the local cell is required to determine the interpolation function of a Pl-Jet. Additional information 
about this interpolation procedure can be found in the original gradient-augmented level set method Q . An 
extension of this nomenclature would be a Full 1-Jet (1-Jet), which would include both the gradient field 
and the first cross-derivatives. Both the Pl-Jet and 1-Jet give similar results and therefore only the Pl-Jet 
is included. Please note that this naming convention differs from Seibold et al. [sj. 

Reinitialization is performed during every time step by explicitly calculating the closest point in a band 
around the interface. This band consists of all nodes within four grid spacings of the interface. A reinitialized 
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Jet is given by 


(l){x) = sgn(.?;>o(a:))| 
'4>{x) = sgn((/)o(a;))- 


X - ccril 
X — Xy 


(44) 

(45) 


||a;-a;r|| 

where (j)[) is the original level set function and xr is the closest point on the interface. The rationale for this 
particular reinitialization scheme over a PDE based one is the direct calculation of the updated gradient 
field. Additionally, it has been shown than the direct reinitialization method is fourth-order accurate for the 
level set field and third-order accurate for the gradient field if the Jet field is smooth If discontinuities 


in the gradient field exist, the accuracy of the level set is locally reduced to second-order while that of the 
gradient to first-order. In a moving interface it is expected that these type of discontinuities are temporally 
temporary. Therefore, this reinitialization should not be detrimental to the overall scheme, which is expected 
to be second-order accurate in space. The use of other schemes, such as that of Cheng and Tsai could 
be extended to reinitialize jet methods, but that is beyond the scope of this work. 

Using Jet schemes there are two choices when calculating curvatures. The first is the direct calculation of 
the curvature at closest points using derivatives of the cubic interpolant. The second is to calculate curvatures 
at grid points and then interpolate to the interface. During the course of the numerical experiments the 
stability of the first method was much worse than the second, and thus the first method was not pursued. 
Therefore, curvatures are first calculated at the grid points by writing Eq. as 


K = 


4^xx4^y 4 “ ^yy4^x ‘^4^xy4^x4^y 


{4>l + 4>l) 


3/2 


(46) 


in two-dimensions with a similar expression for three-dimensions. Eor the 0-Jet all derivatives are computed 
using isotropic finite differences. Eor results using a Pl-Jet the best results were obtained by averaging 
derivatives. Eor example, when computing 4>xx and (jjxy for use in the curvature calculation the following 
averaging would be used: 

<(>. = i , (47) 

4>xx = 5 , (48) 

1 /'n J.x I n „;.M i n 


= I {Dyi)"" -b -b D^y(j)) 


where D^, Dy, D^x and Dxy are the isotropic finite difference approximations for the derivatives and and 
ipy are the components of the gradient field tracked by the Jet. While this does remove some of the locality 
of the Jet method it was found that the accuracy and stability were increased. This also adds additional 
coupling between the components of the Jet. 


4-1. Mean Curvature Flow 

Elow by mean curvature is a standard test problem for assessing the stability and accuracy of interfaces 
moving under a numerically stiff velocity field. In this case the velocity field is given by it = —nn, where n 
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is the mean curvature and n is the outward-facing unit normal vector. 

First consider the benchmark problem of curvature flow applied to a two dimensional circle of radius 
To = 1 centered at (0, 0) in a [—2, 2]^ domain. The initial level set field is given by (f){x, y, 0) = \/x'^ + tq, 
which is a signed distance function. Under mean curvature flow the interface will shrink uniformly where 
the radius at any point in time is given by r{t) = y^rg — 2t. An example of this behavior is provided in Fig. 
|3 (a) I using a grid spacing oi h = 0.125 and a time step of At = 4h^. A comparison of the average radius as 
calculated using the closest points using a first-order and second-order in time scheme is presented in Fig. 
|3(b)[ The second-order method produces more accurate results, as should be expected. 



(a) (b) 

Figure 3: A circle collapsing under mean curvature flow using a Pl-Jet with h = 0.125 and dt = 4h^. (a) The interface location 
is shown at intervals of At until t = 0.375. (b) A comparison between the radius over time using a Pl-Jet with different time 
discretization methods. 


4.1.1. Convergence Study 

The overall accuracy of the method is a combination of the semi-Lagrangian and semi-implicit methods, 

hp-fi 

0{At>^) + 0{^)+0{h^), (50) 

where At is the time step, h is the grid spacing, k is the time step order, p is the order of the interpolant, 
and L is the accuracy order at which the smoothing operator is discretized. Here k = 1 or k = 2, p = 3, and 
L = 2. Clearly, there will be trade-off in accuracy as the time step is reduced. Given a fixed grid spacing 
it should be expected that at large time steps the error will be dominated by the time-integration error, 
0{At^), and as such should decrease at a rate dictated by the time step order. As the time step decreases 
the error begins to be dominated by the smoothing term, 0{h^), it L < p + 1. A further decrease in the 
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time step results in the interpolation error, 0{hP'^^//Si), dominating and any further decrease in the time 
step will result in an increase in the error. 

Therefore there are three convergence regions with respect to the time step which are to be expected. 
The first is a decrease at the rate of O(At^), the second will be a relatively constant error region dictated 
by 0{h^), and the final will be an increase in the error at a rate of Note that the 1/At is 

included to account for the fact that interpolation is performed about 1/At times during the course of a 
simulation. This is not a unique occurrence and has been observed in other instances of semi-Lagrangian 
time integration 
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To verify this accuracy a two-dimensional circle shrinking under mean curvature flow is compared to 
the analytic solution. A circle with an initial radius of 1 centered at the origin is allowed to evolve until 
a time of T = 0.375, at which point the error is calculated. Define the error at a point on the interface 
as e{xT,yT,t) = | -\/a;p -I- j/p — r(t)\, where {xr,yr) is a point on the interface and r{t) = — 2t is the 

analytic radius. The overall error is the Loo-error defined over all closest points calculated at T = 0.375. 
This error definition tests not only the time-evolution of the Jet but also the recovery of the interface by 
way of the cubic interpolation function. 




(a) h = 0.0625 


(b) h = 0.015625 


Figure 4: Comparison of a 0-Jet with a Pl-Jet using both 0{At) and Cl(At^) in time methods for (a) h = 0.0625 and (b) 
h = 0.015625. 


The first error-convergence results are for fixed grid spacings ol h = 0.0625 and h = 0.015625, Fig. |4l 
In this figure both a 0-Jet and a Pl-Jet are used with 0{At) and 0{At^) in time methods. Initially the 
error converges based on the underlying time integration, 0{At^). As the time step decreases the 0{h^) 
terms begin to dominate and approximately constant error regions are observed, particularly for the Pl-Jet 
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using second-order time discretization in Fig. 4(b) As the time step further decreases the dominate term 
becomes 0[h'^/At) and the error increases with a further decrease in the time step. This is clearly seen in 


the error plot of the coarser mesh, Fig. 4(a) Both the 0{At) and 0{At^) methods have roughly the same 
error at small time steps as that error is completely dominated by the interpolation error. For the finer 
mesh, Fig. |4(b)[ the integration error is smaller and therefore this increase at small time steps is delayed. 
In general the PI-Jet provides more accurate results at moderate and large time steps. Only at small time 
steps, when the integration error begins to dominate, does the 0-Jet have slightly smaller error. The exact 
cause of this reversal is not known, but could be due to some fortuitous cancelation. 




(a) 0-Jet 


(b) Pl-Jet 


Figure 5: Error dependence on the time step size for different grid sizes. The time-integration scheme is C)(At^). 


To further explore the influence of the Jet and of the smoothing operator the error as a function of time 
step is investigated for different Jets and grid spacings. Fig. [5] At large time steps the error is completely 
dominated by the time-integration term, and therefore no difference can be observed between the different 
grid spacings, particularly for the Pl-Jet, Fig. |5(b)[ As the grid spacing is decreased the minimum error 
over all time steps decreases at a rate of see Table [T] This indicates that for large time steps any 

spatial errors are dominated by the smoother operator, not the integration error. This is to be expected so 
long as At > h, which results in an overall integration error of 0{h^). 

4-1-2. Cassini Oval 

The extension of the SemiJet method to three dimensions is demonstrated by considering the collapse 
of a Cassini oval under mean-curvature flow, Fig. [SI The initial shape of the interface is given by (/) = 
{{x — a)^ + + a)^ + + z"^) — b'^ with a = 1.29 and b = 1.3 using a grid spacing of h = 0.0634. 
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Table 1: The smallest Loo-error of the interface for the circle under mean curvature flow over all time steps. The overall 
minimum error is dominated by the accuracy of the smoothing operation. 


h 

0-Jet Loo 

0-Jet Order 

Pl-Jet Loo 

Pl-Jet Order 

0.0625 

1.15 X 10-3 

— 

2.43 X 10-3 

—- 

0.03125 

2.45 X 10"^ 

2.23 

5.61 X 10-3 

2.13 

0.015625 

4.03 X 10-3 

2.6 

1.43 X 10-3 

1.96 


This is then replaced by a signed-distance function by reinitialization. The surface of the Cassini moves 
inward and the neck region narrows as the surface shrinks. As time goes on the Cassini pinches off at the 
middle and splits into two separate interfaces. The method has been tested using large time steps and works 
reliably. However, the pinch-off dynamics are very rapid and thus to capture the topological changes the 
time step is chosen to be At = h"^. This example also shows that the use of a PI-Jet does not destroy the 
ability of level set methods to capture topological changes. 

qOq 

(a) t=0 (b) t=0.0605 (c) t=0.0645 (d) t=0.0806 

Figure 6: Collapse of a three-dimensional Cassini oval under mean curvature flow for h = 0.0634 and At = using a Pl-Jet 
with second order time discretization. 




A comparison between the 0-Jet and Pl-Jet for a collapsing Cassini oval is shown in Fig. [71 As can be 
seen the neck using the 0-Jet is thicker than the Pl-Jet at the same time. This indicates that the 0-Jet is 
under-reporting the curvature, leading to a delay in splitting time. 

In Fig. |5]the ability of the Pl-Jet to capture sub-grid information is demonstrated. Here, a comparison is 
made between interfaces constructed with three different interpolation approaches for the collapsing Cassini 
oval at a time of t = 0.0605. The solid black lines represents the interface constructed using the Pl-Jet 
information to determine the cubic interpolant. This result is the same as Fig. |7(b)| This is denoted as 
a sub-grid structure due to the fact that two interface crossings occur on a cell edge where </) > 0 on both 
of the cell corners. Compare this to a cubic interpolant determined using only the level set values and 
approximating the derivatives from finite difference approximations, the dashed line in Fig. [8] In this case 
the center region is completely split. Finally, consider the case of using linear interpolation with only the 
level set values, the dash-dot line in Fig. [5] The description of the interface using linear interpolation is even 
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X 



- 0.2 - 0.1 0 0.1 0.2 


X 


(a) Comparison of a 0-Jet with a Pl-Jet (b) Close-up view of the neck region. 0-Jet under-reports 

the curvature without the higher accuracy afforded by Pl- 
Jet 


Figure 7: Cross-section of the x — y-plane of the Cassini Oval using a 0-Jet and Pl-Jet at a time of f = 0.0605. The neck region 
of the 0-Jet is thicker than the Pl-Jet. 


worse, with only straight lines being described. This demonstrates that the tracking of level set derivatives 
is beneficial in describing the interface. 



Figure 8: Close-up view of the Cassini Oval at t = 0.0605 showing sub-grid structure. The neck region falls within a grid cell 
(represented by solid black lines) which is captured by Pl-Jet through the locally available gradient information 
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4-2. Volume-Conserving Mean Curvature Flow 

In volume-conserving mean curvature flow the velocity of the interface is given by m = — (k — Kavg)'n 
where n is the mean curvature, n is the outward normal to the interface and Kavg is the average mean 
curvature on the interface given by 


/p K da K (5((/)) ||V(/)|| da; 


"^avg 


(51) 


fr IIWII da: 

where fl is the computational domain and S((^) is a Dirac delta function [35j |. Under this velocity field the 
interface will evolve into a single circle in 2D and sphere in 3D, where the average curvature equals the local 
curvature. 






(a) t=0 (b) t=0.10 (c) t=0.12 (d) t=0.15 

Figure 9: Volume-conserving mean curvature flow in two dimensions using a 0-Jet. The grid spacing is h = 0.0634 while the 
time step is At = 0.5h^. At t = 0.15 the volume change is 1.35%. 


The first test case is a two-dimensional system containing two ellipses and a circle subjected to volume 
conserving mean curvature flow. As the curvature is a quantity local to each interface, ideally all interfaces 
would evolve independently until they merge. Provided enough time the final shape of the interface will 
become a circle. Results for both a 0-Jet, Fig. [51 and Pl-Jet, Fig. [151 are provided. The volume loss for 
the 0-Jet is 1.35% while that for the Pl-Jet is 1.28%. 

While the volume conservation is similar it is clear that the 0-Jet result is incorrect. It has been previously 
observed that the non-local nature of the semi-implicit level set method results in the incorrect behavior of 
near-by interfaces 


near each other 


isjl^ 

dm- 


16l| . This is due to the difficulty in calculating the curvature when two interfaces are 


This can be clearly observed in the 0-Jet time-evolution, where the interfaces are 


clearly distorted and have not yet merged at t = 0.15, Fig. [51 Contrast that with the Pl-Jet, where much 
less interface distortion is observed and the interfaces are fully merged at t = 0.15. While the curvature 
calculation is non-local to the interface using the curvature calculation described earlier, the higher accuracy, 
as demonstrated by the results in previous sections, reduces the influence of nearby interfaces. 


A comparison between the original, explicit Jet scheme is also made to the presented SemiJet method by 
considering the volume-conserving mean-curvature flow of a three-dimensional star shape, Fig. 1111 Under 
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(a) t=0 (b) t=0.10 (c) t=0.12 (d) t=0.15 


Figure 10: Volume-conserving mean curvature flow in two dimensions using a PI-Jet. The grid spacing is h = 0.0634 while 
the time step is At = 0.5/i^. At t = 0.15 the volume change is 1.28%. 


this motion any given initial shape will eventually become a sphere. Using the original Jet scheme oscillations 
begin to appear on the interface, which will lead to numerical failure due to incorrect curvature calculation. 
Numerical experiments show that the instability of the original Jet scheme becomes more severe over time 
and the solution eventually breaks down. The presented SemiJet method is able to maintain a smooth 
interface through the entire time with the same time step with only a volume change of 0.37%. 

4-3. Qualitative Stability Analysis 

In this section the stability of the SemiJet method (Pl-Jet with /3 = 0.5) is compared to the original 
semi-implicit level set method (0-Jet with /? = 0.5) and the original Jet scheme (Pl-Jet with /3 = 0). 
Consider a two-dimensional circle under volume-conserving mean curvature flow, u = — (k — Kavg)n. In the 
ideal case the interface will not move due to the local curvature being equal to the average curvature over 
the interface. This situation can be used to investigate the stability of the underlying numerical scheme. 
As the interface should not move any motion of the interface will be due to errors when calculating the 
curvature and advecting the interface. If the time step is too large for a given numerical method these errors 
will compound resulting in the eventual failure of the simulation. 

To test the stability various numerical simulations are carried out to determine the maximum stable 
time step for various combinations of Jet-order, time step accuracy order, and amount of smoothing. While 
the amount of volume change during the course of a simulation is small over reasonable time steps, as 
demonstrated above, during large time steps the volume change can become excessive. To minimize the 
influence of this volume change on the determination of the maximum stable time step a volume correction 
procedure is performed. After every time step the level set is shifted by (14 — Vq) jAn, where 14 is the 
current volume, Vq is the initial volume, and A„ is the current area of the interface. In 2D these are replaced 
with the current area, initial area, and current interface length, respectively. 

A time step is determined to result in a stable simulation if the following conditions are satished: 1) the 
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0 • 


(a) t = 0 


(b) t = 0.05 


(c) t = 0.11 


(d) t = 0.58 



(e)t = 0 (f) 1 = 0.05 (g) 1 = 0.11 (h) 1 = 0.58 

Figure 11: A comparison of the original Jet (/3 = 0) method, (a)-(d), with the Semijet, (e)-(h), showing the volume-conserving 
mean curvature flow of a three dimensional star with Al = 1.5h^ and h = 0.0625. Both solutions use a Pl-Jet. Using the 
original Jet method oscillations on the interface grow over time, leading to the eventual failure of the simulation. The Semijet 
scheme, on the other hand, maintains a smooth solution at all times. Using the Semijet volume loss is about 0.37% at 1 = 0.58. 




simulation can run for 1000 time steps and 2) the difference between the maximum and minimum velocity 
over the last 100 iterations is less than 10% of the minimum velocity over the 1000 total time steps. Sample 
velocity profiles over 1000 iterations for a 0-Jet using /? = 0.5 with second-order accurate time discretization 
for three different time steps is shown in Fig. 1121 There are three types of behaviors observed. The first 
is an unstable time step, where the simulation is not able to complete 1000 iterations before failing due to 
excessive numerical errors. The second type is a semi-stable time step, where the simulation completes 1000 
iterations but the velocity is slowly increasing. This indicates that the simulation will fail at some point in 
the future. Finally, a stable time is one which completes all 1000 iterations and the velocity is stable. A 
stable time step will not have zero velocity; instead any small movement of the interface is adjusted by the 
volume-correction procedure. 

The maximum possible stable time step is reported in Table [2] for a two-dimensional unit circle with 
volume-conserving mean curvature flow on a grid with h = 0.0625. While this determination of a stable 
time step is qualitative in nature it does provide insight into the overall stability properties of the SemiJet 
method. It is clear that the explicit scheme, given by /3 = 0, is restricted by the standard CFL condition of 
At < h'^. The application of a smoothing operator, given by /3 = 0.5, relaxes this constraint. In the case of 
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Figure 12: The maximum velocity using three different time steps using a 0-Jet with /3 = 0.5 for a unit circle under volume- 
conserving mean curvature flow. Three types of behaviors are shown. An unstable time step (At = 1.00) does not complete 
1000 iterations. A semi-stable time step (At = 0.60) completes the 1000 iterations but the velocity is slowly increasing over 
time. A stable time step (At = 0.57) completes 1000 iterations and the velocity remains constant. 


a PI-Jet and first-order time discretization the maximum possible time step is over 400-times that of h?. 

Table 2: The largest stable time step for various combinations of the Jet, time step order, and smoothing on a grid with 
h = 0.0625. The explicit Jet method is given by /3 = 0 while the SemiJet method is given by /3 = 0.5. 



/3 

= 0 

/3 

= 0.5 

Jet Order 

Time Order 

^^max 

I h 

^^max 

^^max ! h 

0 

1 

0.00205 

0.52 

0.98 

251 

PI 

1 

0.00375 

0.96 

1.71 

438 

0 

2 

0.00135 

0.35 

0.57 

146 

PI 

2 

0.00245 

0.63 

0.79 

179 


It is interesting to note that in all cases the use of a Pl-Jet results in a more stable scheme than the 0-Jet. 
This holds true for both the explicit case, ,5 = 0, and the semi-implicit case, /3 = 0.5. It is suspected that the 
averaging of the derivatives when calculating the curvature in the Pl-Jet case is partially responsible for this 
increase in stability. It should also be noted that in general a first-order time discretization scheme is more 
stable than the second-order schemes. This is due to the extrapolation step during the second-order time 
discretization. Any oscillations in the velocity field will be magnified when extrapolating into the future. 
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4-.4- Computational Cost Comparison 

In this section the computational cost of the proposed SemiJet method is compared to a standard fifth 
order upwind WENO method. The first test case considered here is the collapse of a two-dimensional 
unit circle under mean curvature flow, u = —nn. The system is allowed to run until t = 0.375, at which 
time the error is computed, see Sec. 14.11 for more details. The domain spans the [—2,2]^ square. During 
the simulation the wall-clock time spent in the evolution of the level-set or Jet is determined. Note that 
operations common to different cases, such as reinitialization, are not included in the timings. This was 
done to ensure that an accurate comparison between methods is performed. All simulations were performed 
using 8 AMD Opteron 6320 cores. 

The WENO method is discretized using a second-order backward-finite-difference method in time and 
second-order extrapolation of the velocity in time: 


2At 


J- ‘2uji • (pji Uji—i ■ — 0, 


(52) 


where the gradients are discretized using a fifth-order upwind WENO scheme Note that when using 
the WENO scheme, only the level-set field is tracked. The Jet (with /3 = 0) and SemiJet (with /3 = 0.5) 
both utilize a Pl-Jet. As /3 = 0 for the Jet, it is an explicit scheme in time, and thus both the WENO and 
Jet schemes will be under strict stability constraints, see Sec. 14.31 for a qualitative stability analysis of of 
the Jet scheme. It is expected that WENO will obey similar restrictions. The results are presented in Table 
[21 For the WENO and Jet methods the time step is set to 0.25/i^, while with the SemiJet method a larger 
time step can be used. 

It is clear from the results that the Jet and SemiJet are much more costly than the WENO scheme 
for the same grid size and time step, with the Jet being approximately 15 times more expensive and the 
SemiJet approximately 23 times more expensive. This is due to the fact that for every grid point, four sub¬ 
grid points (in 2D) must be updated. Each of these updates requires not only determining the departure 
location, but also the evaluation of multiple cubic interpolants. Note that no effort was made to optimize 
this interpolation, and thus some time savings could be gained there. 

The results also indicate, though, that the quality of the solution, as measured by the Lea error, is much 
better for the Jet and SemiJet schemes as compared to the WENO method. For the same grid size and 
time step, the Jet and SemiJet schemes are between two- and three-orders of magnitude more accurate than 
the WENO scheme. For the SemiJet scheme, this allows for an increase in the time step to still produce an 
accurate result. For example, using a grid size of 257^ (grid spacing of 1.5625 x 10“^) the WENO result has 
an error of 6.81 x 10“^. The same computational cost (^ 100 s) can be achieved with the SemiJet scheme 
using a time step 24 times larger than the WENO case. Even with this larger time step, the error for the 
SemiJet case is three orders of magnitude smaller, 6.82 x 10“® for the SemiJet scheme. Using the SemiJet 
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Table 3: The Loo error and wall-clock time required to model a two-dimensional circle collapsing under mean-curvature flow 
until a time of t = 0.375. 


Method 

Grid Size 

Time Steps 

Isi/h? 

Loo Error 

Elapsed Time (s) 

WENO 

257 

6144 

0.25 

6.81 X 10-3 

99.53 

WENO 

129 

1536 

0.25 

1.33 X 10-2 

6.38 

WENO 

65 

384 

0.25 

2.54 X 10-2 

0.42 

Jet 

257 

6144 

0.25 

2.03 X 10-® 

1527 

Jet 

129 

1536 

0.25 

1.06 X 10-3 

91.95 

Jet 

65 

384 

0.25 

4.68 X 10-3 

5.95 

SemiJet 

257 

6144 

0.25 

2.41 X 10-3 

2377 

SemiJet 

257 

250 

6.14 

6.82 X 10-3 

98.36 

SemiJet 

257 

24 

64 

1.04 X 10-3 

10.50 

SemiJet 

129 

1536 

0.25 

1.94 X 10-3 

148.8 

SemiJet 

129 

48 

8 

2.91 X 10-4 

4.76 

SemiJet 

129 

24 

16 

1.07 X 10-3 

2.45 

SemiJet 

65 

384 

0.25 

6.87 X 10-3 

9.91 

SemiJet 

65 

12 

8 

3.91 X 10-3 

0.35 

SemiJet 

33 

12 

2 

4.33 X 10-3 

0.12 


scheme with a time step 256 times larger (At = 1.563 x 10“^) than the WENO one for the same grid results 
in a scheme which is 6.5 times more accurate but only takes one-tenth of the time, 10.5 s versus 99.53 s. 

This behavior, that at the same computational cost the Semi Jet scheme is much more accurate than the 
WENO scheme is consistent across multiple grid sizes. In fact, using an extremely coarse grid spacing of 
0.125, corresponding to a 33^ grid, with a very large time step still results in a more accurate result than 
the WENO scheme with a fine grid and very small time step, in a computational time which is 830 times 
faster. This behavior is consistent with that shown by Chidyagwai et al. for linear advection Q. 

The second test case is mean curvature flow in three dimensions, u = —0.5/tn, where k represents the 
sum of the principle curvatures. In this case the size of the grid is fixed at 129^ over a cube spanning [—2, 2]^. 
These simulations were performed using 128 AMD Opteron 6320 cores, with the results seen in Table |H 
The results for the three-dimensional case are similar to those for the two-dimensional one. Here, the Jet is 
approximately 26 times more expensive and the SemiJet is approximately 44 times more expensive than the 
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Table 4: The Lao error and wall-clock time required to model a three-dimensional sphere collapsing under mean-curvature 
flow until a time of t = 0.375. All cases use a 129® grid. 


Method 

Time Steps 

At//i2 

Loo Error 

Elapsed Time (s) 

WENO 

1536 

0.25 

2.43 X 10-2 

77.33 

Jet 

1536 

0.25 

1.14 X 10"^ 

2022 

SemiJet 

1536 

0.25 

1.21 X lO-"* 

3417 

SemiJet 

48 

8 

2.73 X 10-4 

116.9 

SemiJet 

30 

12.8 

6.29 X 10-4 

73 

SemiJet 

24 

16 

9.65 X 10-4 

60.15 

SemiJet 

12 

32 

3.63 X 10-3 

32.44 


WENO scheme for the same grid size. In three-dimensions there are eight sub-grid points that need to be 
updated, which explains the increase over the two-dimensional results. Similar to the two-dimensional case, 
the Jet and SemiJet schemes are two orders of magnitude more accurate than the WENO scheme at the 
same grid size and time-step. Using the SemiJet method with 30 time steps, corresponding to At = 0.0125, 
results in the same computational cost as the WENO method but a result which is 40 times more accurate. 
In fact, using an extremely coarse 33^ mesh with a time step of At = 0.125 on two cores results in an elapsed 
time of 21 s and an error of 3.05 x 10“^, which is both faster and more accurate than the WENO case, using 
much fewer computational resources. 

5. Summary and Conclusion 

In this work a robust numerical method, the SemiJet, is presented for moving interface problems with 
numerically-stiff velocity fields. The method builds upon the original Jet scheme developed by Seibold et al. 
and the semi-implicit level set method of Smereka. The influence of the semi-implicit smoothing operation 
is captured by defining a smoothing source term, which is then used when updating the derivative fields in 
the level set Jet. The results demonstrate that the use of a Jet increases both the accuracy and stability of 
the simulations over both the original Jet and original semi-implicit level set methods. Unlike the previous 
SIGALS method only one smoothing operation is required and the use of the e-finite different method allows 
for straight-forward extension to higher order jet schemes. While the SemiJet method is computationally 
more expensive than a WENO scheme for the same grid and time step, the added stability and accuracy 
allows for smaller grid sizes and larger time steps to be used, which results in methods with overall higher 
efficiency and accuracy. 
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